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We show analytically that Newtonian iterations, when applied to 
a polynomial equation, have a positive topological entropy. In a specific ex- 
ample of an attempt to "find" the real solutions of the equation x 1 + 1 = 0, 
we show that the Newton method is chaotic. We analytically find the in- 
variant density and show how this problem relates to that of a piecewise 
linear map. 
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1. Introduction 

Suppose we want to numerically solve a nonlinear equation f(z) = 0, 
where z € R or z £ C. If only calculating the derivative of / is possible, 
perhaps the famous Newton (or Newton-Raphson) method [1] is the first 
method that comes to mind. Starting from some zq, this method uses the 
following iterations: 

Zn+l =Z n ~ = V(Zn) • (1) 

J \ z n) 

If this iteration converges, it usually does so very fast. Because of that, and 
because of its simplicity, the Newton method is one of the most frequently 
used numerical methods ever. Yet it is widely known that even if f(z) is 
a low-order polynomial, boundaries between basins of attractions of different 
roots can be very complicated. Indeed, many pieces of popular software use 
this particular property of the Newtonian iterations to produce aesthetically 
pleasing fractals. 
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Fig. 1. A 2-cycle generated by the Newton method applied to the equation 
x 3 - 2x + 2 = 0. 

The Newton method fails if it hits a zero of the derivative, or when the 
denominator in (pQ) vanishes, unless the zero of the derivative corresponds 
to a multiple root. It is less commonly known that the Newton method 
also fails to find a root if it forms a multicycle. Many textbooks that teach 
the Newton method do not even mention this fact and those that do make 
the impression that multicycles in the Newton method are a rare and unim- 
portant peculiarity. In reality, however, constructing examples that display 
multicycles is quite easy. A reader may verify, for instance, that the points 
{0, 1} form a stable 2-cycle, shown in Fig. [H in the Newton method applied 
to the polynomial x 3 — 2x + 2. Several questions then arise: How many such 
multicycles are there? Is it possible to encounter them in practical numeri- 
cal applications of the Newton method? Is there a connection between the 
multicycles and, say, basins of attraction? 

We may regard consecutive iterations of the Newton method as a dy- 
namical system. Multicycles generated by the Newtonian iterations are the 
periodic orbits of this dynamics. Since the seminal works of Cvitanovic and 
co-workers [2], it is known that the unstable periodic orbits (UPOs) carry 
important information about the dynamical system that generates them. 
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Yet to extract this information, one needs to know at least how many peri- 
odic orbits there are. 

We shall discuss these questions in the case of f(z) being a polynomial. 

This paper is organized as follows: In Section [2] we show that the New- 
tonian iterations applied to a polynomial equation can generate infinitely 
many multicycles and we calculate the topological entropy associated with 
them. In Section [3] we discuss a particularly simple case of the polynomial 
z 2 + l. We show that dynamics resulting from the Newton method, when re- 
stricted to the real axis, is chaotic and equivalent to the dynamics generated 
by a piecewise linear map. We explicitly calculate the invariant density of 
the Newtonian dynamics in this case. Since the number of chaotic systems 
whose invariant densities are known analytically is fairly limited, this result 
is very interesting, at least from a pedagogical point of view. In Section|3]we 
partially generalize these results to the case of polynomials z 2m + 1, m 2. 
Conclusions are given in Section [5j The Appendix contains a proof of the 
theorem that we use to calculate the topological entropy. 



2. Polynomials and multicycles 

Let P(z) be a complex polynomial of order n: 

n 

P(z)=Y,a s z s . (2) 

s=0 

Suppose the polynomial ([2]) has distinct roots^. We solve the polynomial 
equation 

P(z) = (3) 
numerically by Newtonian iterations. This procedure defines a function 

Consecutive iterates of Q satisfy 

^) = ^)-P0^-y * = 1.2,... (5) 

with rp{z) = z. The following Theorem gives a more detailed characteristics 
of the iterates of the function @. It is pivotal in our subsequent discussion: 



Distinct roots of a polynomial must not be confused with roots of muliplicity one. 
For example, the polynomial x 2 (x — l) 2 (a; — 2), of order five, has three distinct roots. 
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Theorem 1. V k ^ 1 

V [Z) P>(z) B k {z) W 

where A k , B k are polynomials of order rejj — rid, where rid is the number of 
distinct roots of the polynomial P(z). 

An elementary, but rather lenghty and technical, proof of this Theorem 
is given in the Appendix. 

Theorem [1] is very useful when we consider multicycles generated by the 
function ([!]). Specifically, all points belonging to a /c-cycle satisfy 

r 1 \z)=z. (7) 

Using Eq. © we get that either P(z) = 0, which means that z is a fixed 
point of the function @, or that 



A k (z) = 0. (8) 

d' "" 1 



This equation has nH — rid solutions on the complex plane. If k is prime, then 




there are N k = (nh — n^)/k different /c-cycles. If k is not prime, we need 
to subtract "spurious" /c-cycles; for example, two rounds over a 2-cycle can 
be mistakenly taken for a 4-cycle. Eventually, for the number of fc-cycles 
that are generated by the Newton method applied to a polynomial with rid 
distinct zeros, we obtain 



"</ i ■ v > • ( 9 ) 

jeD fc / 

where D k is a set of all proper divisors of k (if k is prime, D k = 0). These 
multicycles are periodic orbits of the dynamics generated by successive it- 
erations of the function A limited number of the multicycles may be 
stable, but as in any dynamical system, an overwhelming majority of them 
are unstable and form the unstable periodic orbits (UPOs) of the dynamics. 

The number of multicycles grows very rapidly with k; for example, if 
rid = 5, there are approximately 4.5 • 10 10 17-cycles. To quantify this 
observation, we can calculate the topological entropy [3] of the dynamics 
generated by r/(x) by appropriately counting the UPOs 

S= lim -lniV fe = lnn d . (10) 

fc^oo k 

S > for all rid ^ 2. It is well known that if a map acting on an interval has 
a positive topological entropy, it is chaotic [5]. In a more general setting, 
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if the topological entropy is positive, we can expect some form of chaotic 
behaviour to show up in the dynamics. 

We can see that the Newton method, when applied to a polynomial 
equation, leads to an extremely rich structure of multicycles. But where 
are they? If any of them is stable, as in our example in the Introduction, 
it can show up in a practical application of the Newtonian iterations. The 
unstable majority live on fractal boundaries between basins of attraction. 
It is the dynamics on these boundaries that is the chaotic behaviour whose 
presence is indicated by a positive value of the topological entropy. It is 
also interesting to see that if a polynomial has multiple roots, it has less 
multicycles than a polynomial of the same order and with all roots distinct: 
It is the number of distinct roots, rather than the degree of the polynomial 
or the number of all the roots, that determines the value of the topological 
entropy. 



3. The equation x 2 + 1 = 

In general, multicycles generated by the Newton method live on compli- 
cated geometric objects somewhere on the complex plane. There is, however, 
one nontrivial example where the locale of the multicycles can be pinpointed 
quite accurately. Consider the equation 



z 2 + l = 0. (11) 
It generates the following Newtonian dynamics: 

= 5 (*-;)■ (12) 

It can be proved that if we start with any zq such that Im(zo) > 0, the 
Newtonian iterations converge to = +i. Similarly, if we start with 
Im(zo) < 0, the Newtonian iterations converge to Zoo = —i. Therefore, all 
the multicycles of the function (|12p and all the points that eventually end 
up on them must lie on the real axis. 

We shall, therefore, discuss properties of the map (fT2j) restricted to the 
real axis: 



v ( x ) = - I x - -J , xeR. (13) 

Formally speaking, this map results from an attempt to find the real zeros 
of the equation x 2 + 1 = 0. There are no such zeros, but that does not 
necessarily mean that properties of the map (|13p are not interesting. Fig. [2] 
shows a typical trajectory generated by this map. 
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Fig. 2. A typical trajectory generated by the map (fT3)) . The trajectory spends most 
of its time in the vicinity of x = 0, but from time to time it makes large excursions, 
and then slowly relaxes towards zero. If one waits long enough, arbitrarily large 
excursions can be observed. 



3. 1 . Properties of the map (|13p 
We shall list the most important properties of the map () 1 3 j) : 

• This function is singular in x = and it increases monotonically for 
both x < and x > 0. It does not have a fixed point. 

• Each point has two preimages. The preimages of a point x satisfy 

x ± = x ± \/x 2 + 1 . (14) 

• There are countably many multicycles (UPOs). For each multicycle, 
there are countably many points that eventually fall onto it. 

• x = corresponds to the escape to infinity. There are countably many 
points that are eventually mapped into x = and escape to infinity. 

• The union of all the multicycles, points that lead to them, and points 
that eventually escape to infinity, is also countable and dense in R. 
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This leaves us with uncountably many points that neither escape, nor 
fall on a multicycle, but form a "truly chaotic" orbit of the dynamical 
system (fT3|) . 



Simulations suggest that the map (|13h has an invariant density that 
does not depend on the starting point, provided this point does not escape 
to infinity. This eliminates only countably many points, and in practice 
only three of them: x E { — 1,0,1} (the preimages of ±1 are irrational and 
inaccessible in numerical simulations). Indeed, to compute the invariant 
density, denoted here by p{x), we use the Probenius-Perron equation [BJ: 



p(x)= J 5 (rj(x) - y) p(y) dy . (15) 

— oo 

Using Eq. (|14B . this equation takes the explicit form 



\J x 2 + 1 p(x) = yx + \J x 2 + lj p [x + \J x 2 + 1 

- (x- \Jx 2 + l) p(x- Vx 2 + l) . (16) 

It can be verified by a direct substitution that the Lorentzian distribution 

p{x) = ^ir^y (17) 



solves the Frobenius- Perron equation (|16|) . Thus the map (|13|) joins the 
elite club of maps whose invariant densities are known explicitly. Slowly 
decaying tails of the invariant density (|17p explain the large deviations from 
zero made by the chaotic trajectory. 

3.2. The stereographic projection 

It is interesting to see that the stereographic projection convertss the 
function (|13p into a piecewise linear map on the unit interval. If we substi- 
tute 

x = tan7T0 , (18) 

we obtain 

rj(x) = cot2vr</> = tanvr ( 2(f) + - J . (19) 
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Given the periodicity of the trigonometric function, we can see that the map 
(fT3j) is equivalent to 

C((f>) = 2( t>+\ mod 1 . (20) 

Properties of the map (|20p are easy to find and they nicely illustrate the 
multicycles. 

• The map (|20p is chaotic with the Lyapunov exponent A = In 2. This 
map has a flat invariant density, p((j)) = 1, G [0, 1]. 

• cj) = 1/2 is the fixed point of the map (|20p . The fixed point corresponds 
to an escape to infinity in the language of x given by (|18p . 

• If we represent any <j) G [0, 1] by its binary expansion, the map (|20p 
acts on it by (i) the binary shift, (ii) rejecting the highest bit, and (hi) 
flipping the remaining highest bit. For example, 

0.1100101... 2 -> 0.000101... 2 , (21) 

where the subscript indicates a binary expansion. Thus, all points with 
finite binary expansions eventually converge to O.I2 = 1/2, i.e., escape 
to infinity, points with periodic (starting form some point) binary 
expansions belong to multicycles or settle on one after a finite number 
of steps, and points with nonperiodic, infinite binary expansions lie 
on the chaotic trajectory. Thus the rational numbers either escape 
to infinity under the action of (|20p . or end up on multicycles. As 
the transformation (| 18|) maps rational numbers from the unit interval 
to irrational numbers in M, encountering a true multicycle of (|13p 
during a computer simulation is impossible. However, the fact that 
computer calculations are performed with a finite precision, leads to 
the conclusion that in a computer experiment, the map (|13|) behaves 
as if it were periodic, with a fairly large period, almost regardless of 
the starting point. This is a well known fact about pseudo-random 
sequences. 

• We can count the /c-cycles by counting nontrivially different binary 
sequences of the length k. For example, there are 2 k possibilities of 
distributing {0, 1} among k sites. Two constant sequences must be 
excluded, and as all cyclic permutations are equivalent, we arrive, if k 
is prime, on (2 fc — 2)/k. (If k is not prime, the result of (2 k — 2)/k is 
fractional.) This result, and its generalizations to non-prime fc's, are 
the same as those given by Eq. ([9]) with rid = 2. 
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Fig. 3. The map (fT3"|) (top left) and the map (|20p (bottom left), compared to the 
map (|23| with m = 2 (top right) and the corresponding map on the unit interval, 
obtained numerically (bottom right). 



4. The equation x 2m + 1 = 

The equation 

z 2m + 1 = 0, zeC, m > 2 (22) 

is a natural generalization of the equation (jlip . However, the multicycles 
generated by the Newtonian iterations resulting from the equation (|22p are 
not restricted to the real axis. Indeed, in case of Eq. (|22p the Newton method 
generates identical dynamics on each line e ml on the complex plane, with 
I = 0, 1, . . . m — 1, but there are also fractal boundaries between the basins 
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Fig. 4. The invariant density corresponding to the function (|23p with m = 2, 
obtained numerically. The tails behave like ~ |ir| -4 / 3 , in a full agreement with our 
theoretical predictions. 



of attraction that do not lie on these lines and, apparently, a majority of 
the multicycles are located there. If the Newtonian iterations are restricted 
to the real axis, they are the iterates of the function 

^K^Kw^' x €M. (23) 



This function has qualitatively the same properties as the function (|13[) : it 
is singular at x = 0, grows monotonically from — oo to oo in both domains 
x < and x > 0, and each point has two preimages. Therefore, the number 
of its multicycles (again, on the real axis) must be the same as that of 
(fT3l) . If we use the substitution (fl8l) . the function ([23]) does not convert to 
a piecewise linear map. Fig. [3] compares the maps (fT3j) and ([20]) with the 
map (|23p with m = 2 and the counterpart of C(0)> obtained numerically 
from (g3j| with (ITgjl . 

We now turn to the invariant density. In the asymptotic regime, x>l, 
preimages of a point x satisfy 
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-1 + 

r] (x) ~ x = -x , (24a) 

■q- l (x) ~ x~ = — : -^r tt-. (24b) 

' v ' (2mx) 1 /(2m-i) v ; 

Moreover, if an invariant density exists and is continuous, it must satisfy 

p(x-)~p(0). (24c) 



Thus, if x 3> 1, the Frobenius- Perron equation (|15p takes the approximate 
form 



X2m-1 



2m ( 2mx 

p{x) 



2m — 1 V 2m — 1 



"<°) (25) 



(2m - l)(2m) 1 /(2m-i) ' 



Because of the symmetry of the function (|23|) . the invariant density, if it 
exists, must satisfy p(x) = p(—x). The right-hand side of ([2l)j) is constant. 
Thus, the left-hand side can be constant only if 

p(x) ~ const • | x |- 2 W(2 m -i) ; |a;| _ a, . (26) 

Fig. U] shows the invariant density, found numerically, corresponding to the 
function ()23[) with m = 2. The tails agree perfectly with the theoretical 
prediction (|26|) . Simulations show that for all m ^ 2, the invariant density 
is bimodal, and that the valley between the peaks gets deeper as m increases. 

5. Conclusions 

We have shown that the Newton method, when applied to a polynomial 
equation that has more than one distinct solution, generates infinitely many 
multicycles, and we have calculated the topological entropy associated with 
them. We have also shown that the map (|13p is chaotic, equivalent to 
a piecewise linear map on the unit interval, and we have explicitly found 
its invariant density. This result is interesting by itself as the number of 
maps with analytically known invariant densities is very limited. Because 
the invariant density in question coincides with the Lorentzian probability 
distribution, one could use the function (fT3|) as a basis of a pseudo-random 
generator of the latter. 

From a practical point of view, hitting an unstable multicycle is virtually 
impossible, although it is well known that starting the Newtonian iterations 
in a close vicinity of a fractal boundary between the basins of attractions, 
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where the unstable multicycles live, can slow down the convergence signif- 
icantly. In some cases, stable multicycles can also appear and they can 
prevent the Newton method from converging to a root; the polynomial pre- 
sented in Fig. Q] provides one such example. It would be interesting to see 
what are the criteria for the existence of stable multicycles. In practice, 
however, if a multicycle appears to form, the Newtonian iterations should 
be interrupted by a couple of steps taken with a different method, and with 
the damped Newton method in particular. This usually breaks the multi- 
cycle that has started to form, but can lead to new multicycles in their own 
right. Discussing this point goes beyond the scope of the present paper. We 
should finally mention that the Newton method is not the best method for 
a numerical search for zeros of a polynomial. There are other algorithms, 
better tailored to polynomials. From those the Laguerre method [HE], cou- 
pled with the deflation of the polynomial, is now regarded as the method of 
choice. Despite many arguments in favour of the latter method, the New- 
ton method is probably most commonly used for the task, and this is why 
a study of its UPOs structure is, in our opinion, important. 

To end on a lighter note, we have shown that a numerical search for zeros 
of a polynomial with the Newton method can be quite chaotic, in particular 
when the zeros are just not there. 

LS thanks Mr Wojciech Brzezicki for a helpful discussion. This work 
was supported in part by the Marie Curie Actions Transfer of Knowledge 
project COCOS (contract MTKD-CT-2004-517186). 



Appendix A 

Proof. We give a proof of Theorem Q] from Section [2] in this Appendix. 
First observe that if the polynomial P{z) has multiple roots, there exists 
a polynomial Q(z) such that P(z) = Q(z)R(z) and P'(z) = Q(z)S(z), where 
deg-R(z) = rid, deg5(z) = — 1, is the number of distinct roots of the 
polynomial P(z) and the fraction R(z)/S(z) cannot be further cancelled. 
(If all roots of P{z) are distinct, = n and Q(z) = 1.) Therefore (cf. ©) 

S(z) B k (z) 

To prove the theorem, we need to show that 

I. The polynomials Ak(z), B k {z), generated according to (|5|), have de- 
grees n\ — rid, an d 

II. The polynomials R(z), S(z), A^z), B k (z) do not have common roots, 
from which it follows that the fraction in (jA.ip cannot be further 
cancelled. 
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Both parts of the proof proceed by induction. Note that A\(z) = B\ {£) = 1, 
so the theorem holds for k = 1. For the sake of the notation, let us assume 
that 



rid 



TiZ 3 . 



(A.2) 



Pari I 

Suppose the theorem holds for some k ^ 1 and calculate 



(>)) 



i? z 



R(z) A k {z) 



S(z) B k {z) 



rid / 
r ° + rs ( 



S(z) B k (z) 



(R(z)) 1 



= R(z) 
(note that / - 1 ^ 0) 



1+ e4(0^-')'(^)'^)) 



(S , (*)B fc (*)) n «« 



(S(z)B fc (z)) n " + 



X)r a (5(z)B fc (z)) n -- J] )(-l) l (zS(z)B k (z)r- l (A k (z)) l (R(z)) 1 - 1 

8=1 1 = 1 W 

(A.3) 

Similarly 

R(z) A k {z) 



( z(k) ) 



S z 



S(z) B k (z) 



nd- 1 



(S(z)B k (z))^ 



-Y,( s + Vrs+iiSBkr*- 1 -*^ SB k - RA k ) s . (A.4) 



s=0 



Therefore 



R (zW) _ R{z) N{z) 



S(zW) S(z) B k {z)V(z) 



(A.5) 
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where Af(z), T>(z) are certain polynomials whose forms are implied in 
Eqns. (|A.3|) and (|A.4j) . What are the orders of these polynomials? deg<S = 
n d — 1 and, by assumption, degA k = degB k = — n d . We thus have 

deg(SB k ) n i- s = (n d -s)(n d -l + n k d -n d ) 

= n k d +l -n k d s-n d + s, (A.6) 
deg{zSB k ) s -\A k ) l R 1 - 1 = ( s -l)(l+n d -l+n k -n d ) +l(n k -n d ) + (l-l)n d 

= n k d s-n d . (A.7) 

Thus the term under the sum over s in (|A.3|) is of the order n^ +l — n^s — 
n d + s + n d s = n d +l — 2n + s, which takes a maximal value of n^ +l — n d for 
s = n d . Therefore 

degM(z) =n k+1 -n d . (A.8) 
Similarly, because degzS" = degR = n d , 

degiSB^-'-^zSBk-RAkY = (n d -s-l)(n d -l+n k -n d ) 

+ s(n d +n d -n d ) 
= n k+1 -n k d -n d + s + l, (A.9) 

which takes a maximal value of n^ 1 — for s = n d — 1. Therefore 

degV(z) =n k+1 -n k . (A.10) 

We finally have 

R{r, k (z))_^ R(z)A k (z) R{z) N{z) 



T] k+1 (z) = T] k (z) 



S( V Hz)) S(z) B k {z) S(z) B k (z)V(z) 

R{z) V(z)A k (z) +N(z) _ R{z) A k+1 (z) 
S(z) B k (z)V(z) S(z) B k+l {z) ■ 1 ' ] 

The last statement defines polynomials A k+ i(z), B k+ i{z). Using 
Eqns. (|A,8p and (|A. 1Q|) . it is now easy to verify that degA k+ i = deg-B^+i = 
n^ +1 — n d , provided the last fraction in (jA.lip cannot be cancelled. 
By direct calculations, it is easy to show that 



+ Ep^i^wn^-^^), (A.12.) 

= By(z)S^-Hz)s(z-^^\. (A.12b) 
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Note that we have already shown that A k+ i(z), B k+ \(z) are polynomials. 

Part II 

To complete the proof, we still need to show that the polynomials R(z), 
S(z), A k (z), B k (z) do not have common roots. The polynomials R(z), S(z) 
do not have common roots by construction. Suppose the polynomials R(z), 
S(z), A k (z), B k (z) do not have common roots for some k ^ 1. 

R(z) = 

Let R(z) = 0. Then S(z) / and B k (z) ^ 0, and because of (IA.12bp . 
Bk+i(z) = {B k (z)S{z)) nd / 0. Furthermore, if R{z) = 0, then, by virtue 
of (p£I5a) . A k+ i(z) = A k (z){B k {z)S(z)) n i- l S(z) + 0, because A k (z) + 0, 
either. In other words, R(z) does not have common roots with the other 
polynomials. 

S(z) = 

If S(z) = 0, then z is a zero of the derivative of the original polynomial 
P(z) that is not a multiple root of the latter. In this case the Newton 
iterations diverge and we cannot use the expressions (|A.12p . We need to use 
the definitions impiled in (jA.lip instead; they are still valid by the argument 
of continuity. We have B k +\(z) = Bk(z)T>(z), and if S(z) = 0, this reduces 
to B k +\{z) = s nd ^\(R{z)A k {z)) rid ~ l / 0, where s„ d _i is the highest-order 
coefficient in the polynomial S(z), cf. Eq. (|A.4|) above. For the remaining 
polynomial we have A k +i(z) = A k (z)T>(z) + M(z) which, for S(z) = 0, 
reduces to A k+ i(z) = (-l) nd R nd - 1 (z)Al d (z)(r nd -s nd -i), where r„ d is the 
highest coefficient in R(z). Thus if S(z) = 0, A k+ \{z) could vanish only if 
fn d = s n d -i-, but this is impossible, given the fact that R(z)Q(z) = P(z) 
and S(z)Q(z) = P'(z) for a certain polynomial Q(z). In other words, S(z) 
does not have common roots with the other polynomials. 

Ak+i(z) and Bk+i(z) do not have common roots 

It remains to be shown that A k+ \{z) and B k+ \{z) do not have common 
roots. Let us assume that B k+ \{z) = 0. Because B k+ \(z) = B k (z)V(z), all 
roots of B k (z) are also roots of B k+ \(z). However, if B k (z) = 0, we can 
directly repeat the argument from the preceding paragraph to show that in 
this case A k+1 (z) = {-l) nd R nd -\z)A n k d {z){r nd -s nd ^) + 0. 

Because we have already shown that S(z) does not have common roots 

with B k+ i(z), B k (z) ^ and S (^z — ^{z) ) = i s the only remaining 

possibility for B k+ i(z) to vanish, cf. Eq. (|A.12b|) . If this is the case, then 
form Eq. ()A.12aP we have 
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A k+ i(z) = B^(z)S^{z)R^-\z) R [z - |M . (A.13) 

By assumption, B/,(z) ^ 0, and because we have already shown that 
neither R(z), nor S(z) can have common roots with R(z) ^ 

and S(z) / 0. Finally, because, by construction, the polynomials R(z) and 

S(z) do not have common roots, if S (^z — gj-^fy ) = 0, it follows that 

R ( z ~ W) W)) + °' and therefore ' A k+i{z) + 0. 

We have thus shown that Ak+i(z) and Bk+i(z) do not have common 
roots. This completes the proof. 

□ 
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